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SUMMARY 



A three dimensional finite volume scheme based on hybrid grids containing both tetrahedral and 
hexahedral cells is presented. The application to hybrid grids offers the possibility to combine the 
flexibility of tetrahedral meshes with the accuracy of hexahedral grids. An algorithm to compute a 
dual mesh for the entire computational domain was developed. The dual mesh technique guarantees 
conservation in the whole flow field even at interfaces between hexahedral and tetrahedral domains 
and enables the employment of an accurate upwind flow solver. The hybrid mesh can be adapted 
to the solution by dividing cells in areas of insufficient resolution. The method is tested on different 
viscous and inviscid cases for hypersonic, transsonic and subsonic flows. 


INTRODUCTION 


The development of adaptive methods for efficient and accurate flow calculations has led to two 
major strategies: One class of schemes is based on the application of triangles in two dimensions and 
tetrahedral cells in three dimensions. The employment of those grids for adaptive methods yields some 
advantages: For inviscid calculations the generation of suitable grids even for complex geometries and 
configurations can be done almost automatically by a powerful grid generator. Furthermore, grid 
refinement by introducing new points and retriangulating them can be done by some modules of 
existing grid generators. Such refined grids in general do not require any special treatment either in 
the metrical setup or in the flow solver part of existing codes. Since the grid generation for complex 
configurations has become more and more time consuming in comparison to flow calculation, the first 
point has to be considered to be the main advantage of this class of schemes. 

Besides those typical unstructured methods, other adaptive schemes based on quadrilateral cells, or 
hexahedral cells in three dimensions, exist. Schemes applying to a semi-structured data treatment ([1] 
and [2]) yield only small advantages concerning flexibility when compared with structured methods. 
But even hexahedral schemes with an totally unstructured data treatment, that allow an arbitrary 
cell arrangement, can not compare with tetrahedral schemes in this respect. The main advantage of 
these methods, contrary to tetrahedral schemes, is the higher efficiency and accuracy especially for 
regions dominated by viscosity such as boundary layers in high Reynolds number flows. Cells of high 
aspect ratio can be generated without the need of introducing unwanted small angles. Employing the 
dual mesh technique, control volume faces are orthogonal to the respective edges. 

An approach to combine the advantages of both strategies while circumventing their drawbacks is 
the employment of hybrid grids. The hybrid grids used in the work presented here consist of hexahedral 
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cells in regions near surfaces, where viscous dominated flow can be expected, and of tetrahedral cells 
at some distance from those regions to connect the hexahedral domains and the outer boundaries. 

This paper presents an efficient algorithm to solve the 
Euler- and Navier Stokes equations for hybrid grids con- 
sisting of hexahedral and tetrahedral cells. Because of the 
flexible data structure hexahedral cells do not have to be 
connected to their neighbors in any predefined arrangement 
and almost arbitrary combinations of both cell types are al- 
lowed. In order to obtain a natural coupling between hexahe- 
dral and tetrahedral domains the flow solver part of the pre- 
sented code applies to a dual mesh with control volumes sur- 
rounding each node. Fluxes over the control volume bound- 
aries are determined and related to the respective nodes. 
The dual mesh covers the entire computational domain and 
connects the two grid types. So, the flow solver does not distinguish between the different cell types, 
like some comparable methods ([3]). The use of the dual mesh technique yields automatically conser- 
vation in the flow field and enables the application of an efficient and accurate upwind flow solver to 
calculate the inviscid fluxes. 



node N 0 


METRICAL SETUP 


The input grid data contains information about the geometric positions of the grid nodes and then- 
connections to neighboring nodes. The spatial discretization to determine the viscous and the inviscid 
fluxes for every node is based on the technique of auxiliary cells used as control volumes. The metrical 
setup encloses the determination of size and shape of the dual mesh cells from the given information. 



As shown in figure 1, a control volume surrounding 
node N 0 inside the computational domain has one face 
for every neighboring node 5 . These faces cross the 
edges halfway between N 0 and its neighbors and are 
bounded by the geometric centers of the cells surround- 
ing N 0 . The faces of the control volumes are described 
by face vectors S, representing the face orientation and 
the face size. Each edge of the computational grid has 
one related face in the dual mesh. The face is bounded 
by the midpoints of the cells surrounding the edge. Be- 
cause of the need to deal with hanging nodes, control 
volume faces are not to be bounded by the midpoints of 
Fig. 2: Face of the dual mesh related the faces of computational grid cells as e.g. in [4]. As the 
to an edge connecting N 0 and N, relations between the cells also have to be determined 

in order to connect the midpoints of two neighboring cells, the evaluation of the face vectors for three 
dimensional grids becomes a non trivial task. 
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Consider an edge connecting the nodes N 0 and N x inside the computational domain. This edge is 
shared by n cells C Un as shown in figure 2. The Point N 2 is the midpoint of the edge connecting N 0 
and TVj. The first task in the setup is the determination of the cells C x to C n and the ascertainment 
of their order around the edge. 


The face F is composed of triangles formed by midpoints 
of two neighboring cells and N 2 . The respective face vector 
is the sum of all vectors related to the those triangles. The 
volume f l 0 of the dual cell is composed of tetrahedra with 
the corner points N 0 , N 2 and the midpoints of two neigh- 
boring cells. In a loop running over the cells C x to C n the 
contributions to the face vector and cell volume are evalu- 
ated. 



The metrical setup, including the determination of the 
volume of the dual mesh cells as well as the components 
of the face vectors, has to be executed before the flow cal- 
culation starts. Since the setup forms one dual mesh from 
the entire hybrid grid, conservation is guaranteed for the 

computational domain even at interfaces between regions of hexahedral and tetrahedral cells and at 
interfaces of refined and unrefined cells, as shown in figure 3. 


Fig. 3: Dual mesh at the interface 
between cell types and refinement 
levels 


SPATIAL DISCRETIZATION 
Inviscid Flux Calculation 


The Euler equations are solved employing a modified AUSM Upwind Scheme ([5]) as it is described 
in detail by Kroll and Radespiel ([6]). The special merits of AUSM compared to other upwind schemes 
are the low computational complexity and the low numerical diffusion. 


To compute the inviscid flux over the face F is based on the flow conditions on both sides of the face. 
The values are taken directly from N 0 and N\ for first order calculations. For second order accurate 
calculations the independent flow variables are linearly reconstructed on the control volumes around 
N 0 and N x and 


1 - 

u L = u 0 + Vu 0 ■ -V 01 


( 1 ) 


The gradient Vu 0 of a variable u is obtained by employing a Green-Gauss formula: 


v «o * 7T- ’ £ \ • K + u i) ■ 4 ; ( 2 ) 

;= i z 

where fi 0 is the volume of the dual cell around N 0 and S« is the normal vector of the dual mesh face 
F as shown in figure 1. 
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Near shocks the values on the edges have to be limited to avoid overshoots. The limiting is done 
by an minimum/maximum clipping like it is proposed by Barth in [7j. If a reconstructed value at 
any face of the control volume exceeds the minimum (or maximum) of the values given by node N 0 
and the surrounding nodes N t n , the gradient Vu is scaled by a factor 0, so the reconstructed value 
becomes equal to the minimum (or maximum) of the nodes N 0 n . 


Determination of Viscous Flux 


In order to calculate the viscous flux over the face of an auxiliary cell the primitive variables and 
the derivatives of velocities and temperature on the face have to be computed. For Face F in figure 2 
these values are evaluated by averaging the values of N 0 and N v The gradients of the primitive 
variables in three cartesian directions are the components of Vu. For the calculation of the viscous 
flux the unlimited values have to be used. The determination of the gradients in the way described 
above provides the gradients in the x-, y- and z-direction. This enables the solution of the full Navier 
Stokes equations without relying on a Thin-Layer approximation. 


GRID ADAPTION 


The computational grids can be adapted to the calculated solution by dividing cells in regions of 
insufficient flow resolution. For tetrahedral cells only an isotropic division into eight children cells is 
allowed. Iiexahedral cells can be divided in one, two or all three directions. 



The procedure of tagging cells for division is split into two parts. 
During the first loop the second differences of pressure A p with 


A p — 


Pi ~ 2 Pi + Pz 
Pi + 2 Pi -I- p 3 


Fig. 4: Determination 
helping node N 3 


of 


are computed for each edge. In order to calculate the second dif- 
ferences the most appropriate node N 3 of the neighbors of N t has 
to be determined, as shown in figure 4: N 3 is the node the an- 
gle f3 becomes maximal for. If the value of A p exceeds a certain 
threshold, the cell the respective edge belongs to is tagged for division. For hexahedral cells the rel- 
ative position of the edge in each cell has to be considered in order to divide the cells in the right 
direction. For the second loop all cells that are already tagged are excluded. During this loop other 
flow quantities are computed. In the present code the following criteria are implemented: 
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• first differences of density 

« second differences of density 

• first differences of velocity 


• first differences of velocity weighted by point distance 

If the calculated quantity exceeds a second threshold, the respective cells are also tagged for division. 
The thresholds control the number of divided cells. They can either be fixed before the tagging starts 
or set iteratively in order to receive a predefined number of cells after the cell division is finished. 

For both hexahedral and tetrahedral cells, only one 
hanging node on every edge is allowed. These hanging 
nodes are caused by neighboring divided cells as shown 
in figure 5. The hexahedral cell C t is divided once into 
the children C n and C l2 and introduces a hanging node Fig. 5: Tagging of C 2 in order to 

on the edge of the tetrahedral cell C 2 . If C n was divided P revent a second hanging node 

again, another hanging node on this edge would be introduced. So C 2 has also to be divided before 
the next division of C n . An iterative process runs through the field focussing on the tagged cells. If a 
division would introduce a second hanging node on any edge, the respective cell has also to be tagged 
for division. 



NUMERICAL RESULTS 


The first case is the supersonic inviscid flow (Ma^ = 6) about a blunt body. A coarse hexahedral 
mesh was modified in order to obtain a hybrid mesh with hexahedral cells near the body and tetra- 
hedral cells in some distance from the surface. In this test case the metrical setup and the capability 
of cell division were subject to examinations. Figure 6 shows the position of the shock and the effect 
of a four times refinement that leads to a good shock resolution even for this coarse initial mesh. 

The second case was chosen to test the Navier Stokes formulation of the scheme. For that purpose 
a pure hexahedral 60x40x3 grid over a flat plate was generated. Figure 7 shows the comparison at 
different points between the Blasius solution and the computed solution for a subsonic (Ma^ = 0.5) 
laminar flow with a Reynolds Number of 5000. The grid for this case is not adapted to the solution. 
The results are conform to the analytic solution. Inflow and outflow conditions are computed as 
proposed by Whitfield in [8]. 

The first 3D test case is the transsonic flow around an ONERA M6-Wing. Flow conditions are 
Ma^ = 0.84 with an incidence of 3.06°. The grid has been adapted twice. The initial grid contains 
about 78,000 nodes, 11,000 hexahedra and 370,000 tetrahedra. The once adapted grid contains about 
113,000 nodes, 33,000 hexahedra and 403,000 tetrahedra, while the finest grid contains about 211,000 
nodes, 95,000 hexahedra and 440,000 tetrahedra. The computed flow field is shown in figure 8. The 
characteristic A-shock on the upper side of the wing is nicely resolved. Oscillations at interfaces 
between cells of different refinement levels occure. Those wiggles are subject to investigations in the 
future. 

Also the AGARD 01 test case has been calculated. Figure 9 shows the five times adapted grid and 
the respective solution. The shock resolution is acceptable, the shock regions have been refined during 
each refinement step. The shock on the upper side is located at 0.63 chord length and the shock on 
the lower side at 0.37 chord length. 
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The adaptation for the AGARD 03 test case offers more difficulties. As shown in figure 10 the shock 
behind the airfoil is not resolved very well. The initial grid contains 6,000 nodes. The grid is a three 
times stacked two dimensional grid. The five times adapted grid contains about 55,000 nodes (also 
three stacked planes) and the shock distance is 2.75 chord length behind the trailing edge. 


CONCLUSIONS 


A finite volume scheme using hybrid grids was presented. The employed grids consist of hexahedral 
cells near body surfaces and tetrahedral cells connecting the hexahedral domains and the outer 
boundaries. The use of hexahedral cells offers the possibility to resolve viscous dominated flows such 
as boundary layers efficiently and accurately by applying high aspect ratio cells in those areas. Because 
of the tetrahedral parts, grids become quite flexible and the generation of grids, even for complex 
configurations, is relieved very much compared to structured approaches. 

In a metrical setup a dual mesh is computed from the initial computational grid. This dual mesh 
covers the entire computational domain and connects the two grid types naturally. The feasibility 
of using hybrid grids even for three dimensional flows is shown, but since effective tools for the 
generation of hybrid grids are not available yet at DLR the presented test cases can not prove the 
expected advantages of the approach. The calculation of inviscid fluxes is efficient and accurate. 
Shocks are captured nicely by the employed upwind flow solver. Also the formulation to calculate 
the viscous fluxes has proved its accuracy. Difficulties still occure at interfaces between refined and 
unrefined cells. 

The next step on the way to an automatic system to compute viscous flows around three dimensional 
configurations' is the extension of the hybrid code to prismatic cells in the vicinity of surfaces. The 
prismatic cells can substitute for hexahedral cells, as they yield the same advantages as hexahedral 
cells and algorithms for generating prismatic grids are known from the literature ([9], [10]). The data 
structure of the new version is changed from a point based to an edge based structure. In order to 
enable a vectorization a edge colouring is employed. The performance on a NEC-SX 3 is about 1 
GFlop, including the metrical setup that is not vectorizable. Also the adaption criteria have to be 
improved for the new version. This improvement requires intensive studies of different methods. 


REFERENCES 


[1] M.J. Aftosmis. Viscous flow simulation using an upwind method for hexahedral based adaptive, 
meshes. AIAA 93-0772, 1993. 

[2] J.E. Melton, S.D. Thomas, and G. Cappucccio. Unstructured Euler flow solution using hexahe- 
dral cell refinement. AIAA-91-0637, 1991. 

[3] M. Soetrisno, S.T. Imlay, and D.W. Roberts. A zonal implicit procedure for hybrid structured- 
unstructured grids. AIAA-94-0645, 1994. 


176 



[4] P. Vijayan and Y. Kallinderis. A 3d finite volume scheme for the Euler equations on adaptive 
tetrahedral grids. Journal of Computational Physics, 113:249-267, 1994. 

[5] M.S. Liou and C.J. Steffen. A new flux splitting scheme. Computers and Fluids, 107(l):23-39, 
1993. 

[6] N. Kroll and R. Radespiel. An improved flux vector split discretisation scheme for viscous flows. 
DLR-FB 93-53, 1993. 

[7] T.J. Barth and D.C. Jesperson. The design and application of upwind schemes on unstructured 
meshes. AIAA-89-0366, 1989. 

[8] D.L. Whitfield. Three-dimensional unsteady Euler equation solutions using flux vector splitting, 
1983. 

[9] K. Nakahashi. Adaptive prismatic grid method for external viscous flow computations. AIAA- 
93-3314-CP, 1993. 

[10] M.J. Marchant and N.P. Weatherill. Unstructured grid generation for viscous flow simulations. In 
N.P. Weatherill, P.R. Eiseman, J. Hauser, and J.F. Thompson, editors, Numerical grid generation 
in computational fluid dynamics and related fields, pages 151-162, 1994. 



Fig. 6: Flow around a blunt body: Initial grid and four times adapted grid with respec- 
tive solutions 
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Fig. 7: Laminar flow over a flat plate: Computed Solution compared to analytic Solution 
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Fig. 9: AGARD01 test case: Five times adapted grid and respective solution 



Fig. 10: AGARD03 test case: Initial grid and initial solution and five times adapted 
grid with respective solution 



